library(MASS)

setwd("C:\\Users\\mdw\\Desktop\\social distance september 2007\\Data")


dat<-dget("bcow_052007")
dattxt<-dget("bcowtxt_052007")

kncaus<-subset(dat, dat$bosnia==0)
kncaustxt<-subset(dattxt, dattxt$bosnia==0)

kbos<-subset(dat, dat$bosnia==1)
kbostxt<-subset(dattxt, dattxt$bosnia==1)

choose.vars<-c("religiosity", 
               "nat.feeling", 
               "representation", 
               "ethnicfriends", 
               "trust",           
               "history",           
               "nat.trust",         
               "pride",            
               "newcomers",         
               "separation")
        
knauchose<-kncaus;  knauchose$"nationality"<-kncaustxt$"nationality"
kboschose<-kbos;    kboschose$"nationality"<-kbostxt$"nationality"

slect.bos<-subset(kboschose, select=c("nationality", "region2", choose.vars))
slect.nca<-subset(knauchose, select=c("nationality", "region2", choose.vars))

slect.bos$nationality<-as.character(slect.bos$nationality)
slect.nca$nationality<-as.character(slect.nca$nationality)

slect.bos[slect.bos==9]<-NA
slect.bos[slect.bos==8]<-NA
slect.nca[slect.nca==9]<-NA
slect.nca[slect.nca==8]<-NA

slect.bos<-na.omit(slect.bos) 
slect.nca<-na.omit(slect.nca) 

## plots by ethnicity:

# serbs = red
# croats = blue
# bosnian = green 
# the rest = gray

color.bos<-slect.bos$nationality
color.bos[color.bos!="bosniac"&color.bos!="serb"&color.bos!="croat"]<-"gray"
color.bos[color.bos=="serb"]<-"red"
color.bos[color.bos=="croat"]<-"blue"
color.bos[color.bos=="bosniac"]<-"green"

# Russians = red
# Ossetins = green
# avar = blue
# the rest = gray

color.nca<-slect.nca$nationality
color.nca[color.nca!="russian"&color.nca!="avar"&color.nca!="ossetin"&color.nca!="kabardinian"&color.nca!="darghinn"]<-"gray"

color.nca[color.nca=="russian"]<-"red"
color.nca[color.nca=="avar"]<-"blue"
color.nca[color.nca=="ossetin"]<-"green"
color.nca[color.nca=="kabardinian"]<-"orange"
color.nca[color.nca=="darghinn"]<-"purple"


### 

# Bosnia first: 
d<-as.matrix(dist(slect.bos[,-c(1,2)], diag=T))
par(mfrow=c(1,1))
loc <- cmdscale(d)
     x <- -loc[,1]
     y <- loc[,2]
     plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, col=color.bos, xaxt="n",yaxt="n")#,bty="n", xlim=c(-10, 9), ylim=c(-8,8))#; abline(v=c(-2, 9))
leg.txt<-c("Bosniac ", "Serbian ", "Croat ", "All others")
legend(locator(n=1), leg.txt, pch =NULL, bty = "n",
       text.col=c("green", "red", "blue", "gray"), cex = 1.1)


### Caucasus
d.nca<-as.matrix(dist(slect.nca[,-c(1,2)], diag=T))
par(mfrow=c(1,1))
loc <- cmdscale(d.nca, k=3)
     x <- loc[,1]
     y <- loc[,2]

par(mfrow=c(1,1))
     plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, 
          col=color.nca, xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)),
          ylim=c(min(y),max(y)));# abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=3, lty=2)
leg.txt<-c("Russian", "Avar", "Ossetin", "Kabardin", "Dargin", "All others")
legend(locator(n=1), leg.txt, pch =NULL, bty = "n",
       text.col=c("red", "blue", "green", "orange", "purple", "gray"), cex = 1.1)
       
       
# color scheme for regions:
color.bos<-slect.bos$region2
color.bos[color.bos==1]<-"green"
color.bos[color.bos==2]<-"red"
color.bos[color.bos==3]<-"blue"

color.nca<-slect.nca$region2
color.nca[color.nca==42]<-"red"
color.nca[color.nca==43]<-"orange"
color.nca[color.nca==46]<-"blue"
color.nca[color.nca==47]<-"gray"
color.nca[color.nca==45]<-"green"

# Bosnia by region: 
d<-as.matrix(dist(slect.bos[,-c(1,2)], diag=T))
par(mfrow=c(1,1))
loc <- cmdscale(d)
     x <- loc[,1]
     y <- loc[,2]
     plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, col=color.bos, xaxt="n",yaxt="n", xlim=c(min(x)-1, max(x)+1), ylim=c(min(y),max(y+2)))#; abline(v=c(-2, 9))
leg.txt<-c("Federation of Bosnia and Herzegovina", "Republika Srpska", "Brcko Federal District")
legend(locator(n=1), leg.txt, pch =NULL, bty = "n",
       text.col=c("green", "red", "blue"), cex = 1.1)

# Caucasus by region:
d.nca<-as.matrix(dist(slect.nca[,-c(1,2)], diag=T))
par(mfrow=c(1,1))
loc <- cmdscale(d.nca, k=3)
     x <- loc[,1]
     y <- loc[,2]
     plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, col=color.nca, xaxt="n",yaxt="n", xlim=c(min(x)-1, max(x)+1), ylim=c(min(y),max(y+2)))#; abline(v=c(-2, 9))
leg.txt<-c("Stavropol", "Karachay Cherkessia", "Dagestan", "Kabardino Balkaria", "North Ossetia")
legend(locator(n=1), leg.txt, pch =NULL, bty = "n",
       text.col=c("red", "orange", "blue", "gray", "green"), cex = 1.3)
      

## need to fix the color for the kabardinian and darghinn
par(mfrow=c(2,3))
     plot(x[slect.nca$nationality=="russian"], 
          y[slect.nca$nationality=="russian"], 
          xlab="", ylab="", main="russian", pch=19, bty="n", cex=1.5, col="red", 
          xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)));
          abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)
          
     plot(x[slect.nca$nationality=="ossetin"], 
          y[slect.nca$nationality=="ossetin"], 
          xlab="", ylab="", main="ossetin", pch=19, bty="n", cex=1.5, col="green", 
          xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
          abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)
          
     plot(x[slect.nca$nationality=="avar"], y[slect.nca$nationality=="avar"], 
          xlab="", ylab="", main="avar", pch=19, bty="n", cex=1.5, col="blue", 
          xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
          abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)
          
     plot(x[slect.nca$nationality=="kabardinian"], y[slect.nca$nationality=="kabardinian"], 
          xlab="", ylab="", main="kabardin", pch=19, bty="n", cex=1.5, col="orange", 
          xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
          abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)
          
     plot(x[slect.nca$nationality=="darghinn"], y[slect.nca$nationality=="darghinn"], 
          xlab="", ylab="", main="dargin", pch=19, bty="n", cex=1.5, col="purple", 
          xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
          abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)
          
     plot(x[slect.nca$nationality!="russian"&
            slect.nca$nationality!="ossetin"&
            slect.nca$nationality!="avar"&slect.nca$nationality!="kabardinian"&slect.nca$nationality!="darghinn"], 
            y[slect.nca$nationality!="russian"&
            slect.nca$nationality!="ossetin"&
            slect.nca$nationality!="avar"&slect.nca$nationality!="kabardinian"&slect.nca$nationality!="darghinn"], 
            xlab="", ylab="", main="rest", pch=19, bty="n", cex=1.5, col="gray", 
            xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
            abline(v=c(mean(x)), h=c(mean(y)), col="purple", lwd=2, lty=2)



##    color latent positions by questions answered: a loop to do postscripts
##    the color scheme is quite arbitrary. do N Caucasus first. 
d.nca<-as.matrix(dist(slect.nca[,-c(1,2)], diag=T))
loc <- cmdscale(d.nca, k=3)
     x <- loc[,1]
     y <- loc[,2]

Save<-NULL
for (i in 1:(dim(slect.nca)[2]-2)) {
color.q<-slect.nca[, 2+i]
scal<-sort(unique(slect.nca[, 2+i]))

if (length(scal)==2){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==3){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-1]<-"lightgreen"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==4){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-2]<-"green"
                     color.q[color.q==max(scal)-1]<-"gray"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==5){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-3]<-"green"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==6){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-4]<-"lightblue";
                     color.q[color.q==max(scal)-3]<-"green"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==7){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-5]<-"lightblue";
                     color.q[color.q==max(scal)-4]<-"green"
                     color.q[color.q==max(scal)-3]<-"lightgreen"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
postscript(paste(colnames(slect.nca)[i+2], ".ps", sep=""), horizontal=F,height=7,width=7)
par(mfrow=c(1,1))
plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, 
col=color.q, xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
dev.off()
Save<-c(Save, paste(colnames(slect.nca)[i+2], ".ps", sep=""))
}


d<-as.matrix(dist(slect.bos[,-c(1,2)], diag=T))
loc <- cmdscale(d)
     x <- loc[,1]
     y <- loc[,2]
         
Save<-NULL
for (i in 1:(dim(slect.bos)[2]-2)) {
color.q<-slect.bos[, 2+i]
scal<-sort(unique(slect.bos[, 2+i]))

if (length(scal)==2){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==3){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-1]<-"lightgreen"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==4){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-2]<-"green"
                     color.q[color.q==max(scal)-1]<-"gray"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==5){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-3]<-"green"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==6){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-4]<-"lightblue";
                     color.q[color.q==max(scal)-3]<-"green"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
if (length(scal)==7){color.q[color.q==min(scal)]<-"blue";
                     color.q[color.q==max(scal)-5]<-"lightblue";
                     color.q[color.q==max(scal)-4]<-"green"
                     color.q[color.q==max(scal)-3]<-"lightgreen"
                     color.q[color.q==max(scal)-2]<-"gray"
                     color.q[color.q==max(scal)-1]<-"orange"
                     color.q[color.q==max(scal)]<-"red"}
                     
postscript(paste(colnames(slect.bos)[2+i], "_b.ps", sep=""), horizontal=F,height=7,width=7)
par(mfrow=c(1,1))
plot(x, y, xlab="", ylab="", main="", pch=19, bty="n", cex=1.5, 
col=color.q, xaxt="n",yaxt="n",bty="n", xlim=c(min(x), max(x)), ylim=c(min(y),max(y)))
dev.off()
Save<-c(Save, paste(colnames(slect.bos)[i+2], ".ps", sep=""))
}


### the graph/figure part ends here ::::


## calculate mean and dv of social distances:
 
# for ethnicity:
mean(d[slect.bos$nationality=="croat", slect.bos$nationality=="croat"])  
sd(as.vector(d[slect.bos$nationality=="croat", slect.bos$nationality=="croat"]))

mean(d.nca[slect.nca$nationality=="kabardinian", slect.nca$nationality=="kabardinian"])
mean(d.nca[slect.nca$nationality=="russian", slect.nca$nationality=="avar"])
mean(d.nca[slect.nca$nationality=="russian", slect.nca$nationality=="ossetin"])
mean(d.nca[slect.nca$nationality=="avar", slect.nca$nationality=="ossetin"])


# for regions:
mean(d[slect.bos$region2==1, slect.bos$region2==1])  
sd(as.vector(d[slect.bos$region2==1, slect.bos$region2==1]))
mean(d[slect.bos$region2==2, slect.bos$region2==2])  
sd(as.vector(d[slect.bos$region2==2, slect.bos$region2==2]))
mean(d[slect.bos$region2==3, slect.bos$region2==3])  
sd(as.vector(d[slect.bos$region2==3, slect.bos$region2==3]))

mean(d[slect.bos$region2==1, slect.bos$region2==2])  
mean(d[slect.bos$region2==1, slect.bos$region2==3])  
mean(d[slect.bos$region2==2, slect.bos$region2==3])  

dim(subset(slect.bos, slect.bos$region2==2))
dim(subset(slect.nca, slect.nca$region2==42))


mean(d[slect.nca$region2==42, slect.nca$region2==42])  
sd(as.vector(d[slect.nca$region2==42, slect.nca$region2==42]))

mean(d[slect.nca$region2==43, slect.nca$region2==43])  
sd(as.vector(d[slect.nca$region2==43, slect.nca$region2==43]))

mean(d[slect.nca$region2==45, slect.nca$region2==45])  
sd(as.vector(d[slect.nca$region2==45, slect.nca$region2==45]))

mean(d[slect.nca$region2==46, slect.nca$region2==46])  
sd(as.vector(d[slect.nca$region2==46, slect.nca$region2==46]))

mean(d[slect.nca$region2==47, slect.nca$region2==47])  
sd(as.vector(d[slect.nca$region2==47, slect.nca$region2==47]))


mean(d[slect.nca$region2==42, slect.nca$region2==43])  

mean(d[slect.nca$region2==42, slect.nca$region2==45]) 

mean(d[slect.nca$region2==42, slect.nca$region2==46])  

mean(d[slect.nca$region2==42, slect.nca$region2==47])  

mean(d[slect.nca$region2==43, slect.nca$region2==45])  

mean(d[slect.nca$region2==43, slect.nca$region2==46])  

mean(d[slect.nca$region2==43, slect.nca$region2==47])  

mean(d[slect.nca$region2==45, slect.nca$region2==46])  

mean(d[slect.nca$region2==45, slect.nca$region2==47])

mean(d[slect.nca$region2==46, slect.nca$region2==47])



## region names for North Caucasus:

Federation of Bosnia and Herzegovina                     Republika Srpska 
                                1114                                  802 
              Brcko Federal District                            Stavropol 
                                  84                                  810 
                 Karachay Cherkessia                             Dagestan 
                                 121                                  625 
                  Kabardino Balkaria                        North Ossetia 
                                 246                                  198 
table(dat$region2)

   1    2    3   42   43   45   46   47 
1114  802   84  810  121  625  246  198 
